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Abstract 

In this letter, we note that the denoising performance of Non-Local 
Means (NLM) at large noise levels can be improved by replacing the 
mean by the Euclidean median. We call this new denoising algorithm 
the Non-Local Euclidean Medians (NLEM). At the heart of NLEM 
is the observation that the median is more robust to outliers than 
the mean. In particular, we provide a simple geometric insight that 
explains why NLEM performs better than NLM in the vicinity of 
edges, particularly at large noise levels. NLEM can be efficiently im- 
plemented using iteratively reweighted least squares, and its compu- 
tational complexity is comparable to that of NLM. We provide some 
preliminary results to study the proposed algorithm and to compare 
it with NLM. 

Keywords: Image denoising, non-local means, Euclidean median, iter- 
atively reweighted least squares (IRLS), Weiszfeld algorithm. 

1 Introduction 

Non-Local Means (NLM) is a data-driven diffusion mechanism that was 
introduced by Buades et al. in HI. It has proved to be a simple yet powerful 
method for image denoising. In this method, a given pixel is denoised 
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using a weighted average of other pixels in the (noisy) image. In particular, 
given a noisy image u = {ui), the denoised image u = {ui) at pixel i is 
computed using the formula 

Ui = ^ , (1) 

where Wij is some weight (affinity) assigned to pixels i and j. The sum 
in ^ is ideally performed over the whole image. In practice, however, 
one restricts j to a geometric neighborhood of i, e.g., to a sufficiently large 
window of size S x S 

The central idea in ||T| was to set the weights using image patches cen- 
tered around each pixel, namely as 

Wij =exp(^- ^\\Pi-Pjfy (2) 

where Pj and Pj are the image patches of size k x k centered at pixels i 
and j. Here, ||P|| is the Euclidean norm of patch P as a point in R'^ , and 
/i is a smoothing parameter. Thus, pixels with similar neighborhoods are 
given larger weights compared to pixels whose neighborhoods look differ- 
ent. The algorithm makes explicit use of the fact that repetitive patterns 
appear in most natural images. It is remarkable that the simple formula 
in ^ often provides state-of-the-art results in image denoising [2J. One 
outstanding feature of NLM is that, in comparison to other denoising tech- 
niques such as Gaussian smoothing, anisotropic diffusion, total variation 
denoising, and wavelet regularization, the so-called method noise (differ- 
ence of the denoised and the noisy image) in NLM appears more like white 
noise HHH. We refer the reader to [2l for a detailed review of the algorithm. 

The rest of the letter is organized as follows. In Section |2l we explain 
how the denoising performance of NLM can be improved in the vicinity of 
edges using the Euclidean median. Based on this observation, we propose a 
new denoising algorithm in Section|3]and discuss its implementation. This 
is followed by some preliminary denoising results in Section HI where we 
compare our new algorithm with NLM. We conclude with some remarks 
in Section |5l 



2 Better robustness using Euclidean median 

The denoising performance of NLM depends entirely on the reliability of 
the weights Wij . The weights are, however, computed from the noisy image 
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and not the clean image. Noise affects the distribution of weights, particu- 
larly when the noise is large. By noise, we will always mean zero-centered 
white Gaussian noise with variance in the rest of the discussion. 

To understand the effect of noise and to motivate the main idea of the 
paper, we perform a simple experiment. We consider the particular case 
where the pixel of interest is close to an (ideal) edge. For this, we take 
a clean edge profile of unit height and add noise {a = 0.2) to it. This is 
shown in Fig. [TJ We now select a point of interest a few pixels to the right of 
the edge (marked with a star). The goal is to estimate its true value from its 
neighboring points using NLM. To do so, we take 3-sample patches around 
each point {k = 3), and a search window of 5 = 41. The patches are shown 
as points in 3-dimensions in Fig|2al The clean patch for the point of interest 
is at (1,1,1). 

We now use (|2]) to compute the weights, where we set h = Wa. The 
weights corresponding to the point of interest are shown in Fig. [ibl Using 
the noisy weights, we obtain an estimate of around 0.65. This estimate 
has a geometric interpretation. It is the center coordinate of the Euclidean 
mean WjPj/ Wj, where wj are the weights in Fig. [ibl and are the 
patches in Fig. |2al The Euclidean mean is marked with a purple diamond in 
Fig. 121 Note that the patches drawn from the search window are clustered 
around the centers A = (0, 0, 0) and B = (1, 1, 1). For the point of interest, 
the points around A are the outliers, while the ones around B are the inliers. 
The noise reduces the relative gap in the corresponding weights, causing 
the mean to drift away from the inliers toward the outliers. 

Note that the Euclidean mean is the minimizer of Yj Wj\\P — PjW^ over 
all patches P. Our main observation is that, if we instead compute the 
minimizer of J2j ~ Pjll over all P, and take the center coordinate of 

P, then we get a much better estimate. Indeed, the denoised value turns out 
to be around 0.92 in this case. The above minimizer is called the Euclidean 
median (or the geometric median) in the literature O. We will often simply 
call it the median. We repeated the above experiment using several noise 
realizations, and consistently got better results using the median. Averaged 
over 10 trials, the denoised value using the mean and median were found 
to be 0.62 and 0.93, respectively. Indeed, in Fig. |2l note how close the 
median is to the true patch compared to the mean. This does not come as 
a surprise since it is well-known that the median is more robust to outliers 
than the mean. This fact has a rather simple explanation in one dimension. 
In higher dimensions, note that the former is minimizer of (the square of) 
the weighted P norm of the distances j|P — Pj\\, while the latter is the 
minimizer of the weighted norm of these distances. It is this use of the 
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(b) Weights. 
Figure 1: Ideal edge in one dimension. 
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(b) 2d projection (first 2 coordinates). 

Figure 2: Outlier model of the patch space for the point of interest in Fig. 
Hal Due to its robustness to outliers, the Euclidean median behaves as a 
better estimator of the clean patch than the Euclidean mean (see text for 
explanation). 
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norm over the i'^ norm that makes the Euclidean median more robust to 
outliers ||3l. 



3 Non-Local Euclidean Medians 

Following the above observation, we propose Algorithm[T]below which we 
call the Non-Local Euclidean Medians (NLEM). We use S{i) to denote the 
search window of size 5x5 centered at pixel i. 

Algorithm 1 Non-Local Euclidean Medians 

Input: Noisy image u = (ui), and parameters h, A, S, k. 
Return: Denoised image u = (ui). 

(1) Estimate noise variance cr^, and set h = Xa 

(2) Extract patch Pj G R*^ at every pixel i. 

(3) For every pixel i, do 

(a) Set Wij = exp(-||Pi - Pj|p//i^) for every j G S{i). 

(b) Find patch P that minimizes J2jeS{i) 11-^ ~ II- 

(c) Assign Ui the value of the center pixel in P. 



The difference with NLM is in step (3)b which involves the computa- 
tion of the Euclidean median. That is, given points xi, . . . ,Xn € R'^ and 
weights wi, . . . , Wn, we are required to find x G R'^ that minimizes the con- 
vex cost X^"=i Wj\\x — Xj\\. There exists an extensive literature on the com- 
putation of the Euclidean median; see 0113/ and the references therein. The 
simplest algorithm in this area is the so-called Weiszfeld algorithm |4l|5l. 
This is, in fact, based on the method of iteratively reweighted least squares 
(IRLS), which has received renewed interest in the compressed sensing 
community in the context of minimization [61 [Z|- Starting from an es- 
timate x^^\ the idea is to set the next iterate as 



x^ ' = are; mm > Wi- — ttt — - — -. 

cceR''^ ^ kE('=) - xA\ 
j=i II ^11 

This is a least-squares problem, and the minimizer is given by 

,(fc+i) _ ^J^i 



= (3) 



where ^- = Wj/\\x'^'''> — Xj\\. Starting with an initial guess, one keeps 
repeating this process until convergence. In practice, one needs to address 
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the situation when x^''^ gets close to some Xj, which causes fij to blow 
up. In the Weiszfeld algorithm, one keeps track of the proximity of x^^^ to 
all the Xj, and a;'^'^"''^) is set to be Xi if \\x^'''^ — Xi\\ < e for some i. It has 
been proved by many authors that the iterates converge globally, and even 
linearly (exponentially fast) under certain conditions, e.g., see discussion in 

Following the recent ideas in [Zl, we have also tried regularizing Q 
by adding a small bias, namely by setting /ij'^'' = Wj/{\\x'^^^ — XjW^ + ^V)^^"^- 
The bias eu is updated at each iterate, e.g., one starts with eo = 1 and grad- 
ually shrinks it to zero. The convergence properties of a particular flavor 
of this algorithm are discussed in |l7|. We have tried both the Weiszfeld al- 
gorithm and the one in 161. Experiments show us that faster convergence 
is obtained using the latter, typically within 3 to 4 steps. The overall com- 
putational complexity of NLEM is k"^ ■ S"^ ■ Nnev per pixel, where A'^itcr is the 
average number of IRLS iterations. The complexity of the standard NLM is 
k'^ ■ S'^ per pixel. 

4 Experiments 

We now present the result of some preliminary denoising experiments. For 
this, we have used the synthetic images shown in Fig. |3] and some nat- 
ural images. For NLEM, we computed the Euclidean median using the 
simple IRLS scheme described in [6J. For all experiments, we have set 

5 = 21, k = 7, and h = Wa for both algorithms. These are the standard 
settings originally proposed in [IJ. 

First, we considered the Checker image. We added noise with a = WO 
resulting in a peak-signal-to-noise ratio (PSNR) of 8.18dB. The PSNR im- 
proved to 17.94dB after applying NLM, and with NLEM this further went 
up to 19.45dB (averaged over 10 noise realizations). This 1.5dB improve- 
ment is perfectly explained by the arguments provided in Section|2l Indeed, 
in Fig nil we have marked those pixels where the estimate from NLEM is 
significantly closer to the clean image than that obtained using NLM. More 
precisely, denote the original image by fi, the noisy image by Uj, and the 
denoised images from NLM and NLEM by Ui and u'^. Then the "+" in the 
figure denotes pixels i where |n- — | < \ui — fi\ — 10. Note that these points 
are concentrated around the edges where the median performs better than 
the mean. 

So what happens when we change the noise level? We went from cr = 10 
to cj = 100 in steps of 10. The plots of the corresponding PSNRs are shown 
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(a) Checker ((t=100). 



(b) House (cr=60). 



Figure 4: The "+" symbol marks pixels where the estimate returned by 
NLEM is significantly better than that returned by NLM. 
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in Fig. |5al At low noise levels (a < 30), we see that NLM performs as good 
or even better than NLEM. This is because at low noise levels the true neigh- 
bors in patch space are well identified, at least in the smooth regions. The 
difference between them is then mainly due to noise, and since the noise 
is Gaussian, the least-squares approach in NLM gives statistically optimal 
results in the smooth regions. On the other hand, at low noise levels, the 
two clusters in Fig. |2]are well separated and hence the weights Wij for NLM 
are good enough to push the mean towards the right cluster. The median 
and mean in this case are thus very close in the vicinity of edges. At higher 
noise levels, the situation completely reverses and NLEM performs consis- 
tently better than NLM. In Fig. |5al we see this phase transition occurs at 
around a = 30. The improvement in PSNR is quite significant beyond this 
noise level, ranging between 0.5dB and 2.1dB. The exact PSNRs are given 
in Table [H 

Next, we tried the above experiment on the Circles image. The results 
are given in table [H The phase transition in thsi case occurs around a = 25. 
NLM performs significantly better before this point, but beyond the phase 
transition, NLEM begins to perform better, and the gain in PSNR over NLM 
can be as large as 2.2dB. The method noise for NLM and NLEM obtained 
from a typical denoising experiment are shown in Fig. |6l Note that the 
method noise for NLEM appears more white (with less structural features) 
than that for NLM. 

Finally, we considered some benchmark natural images, namely House, 
Barbara, and Lena. The PSNRs obtained from NLM and NLEM for these 
images at different noise levels are shown in Table [H The table also com- 
pares the Structural SIMilarity (SSIM) indices |1T2| obtained from the two 
methods. Note that a phase transition in SSIM is also observed for some of 
the images. In Fig. |4R we show the pixels where NLEM does better (in the 
sense defined earlier) than NLM for the House image at a = 60. We again 
see that these pixels are concentrated around the edges. 

The improvements in PSNR and SSIM are quite dramatic for the syn- 
thetic images Checker and Circles. This is expected because they contain 
many edges and the edges have large transitions. The improvement in 
PSNR and SSIM are less dramatic for the natural images. But, consider- 
ing that NLM already provides top quality denoising results, this small 
improvement is already significant. We have also noticed that the perfor- 
mance of NLEM (and that of NLM) can be further improved by considering 
only the top 50% of the 5^ neighbors with largest weights. This simple mod- 
ification improves both NLM and NLEM (and also their run time), while 
still maintaining their order in terms of PSNR performance. A compari- 
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Figure 5: Comparison of denoising performance for the Checker image at 
noise a = 80. NLM-kNN and NLEM-kNN refer to the respective modifica- 
tions of NLM and NLEM where we only use the top 50% of the weights. 
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Figure 6: Comparison of the method noise for the Circles image at cr = 80. 

son of the four different methods is given in Fig. |5bl The results provided 
here are rather incomplete, but they already provide a good indication of 
the denoising potential of NLEM at large noise levels. While we have con- 
sidered a fixed setting for the parameters, our general observation based 
on extensive experiments is that NLEM consistently performs better than 
NLM beyond the phase transition, irrespective of the parameter settings. In 
future, we plan to investigate ways of further improving NLEM, and study 
the effect of the parameters on its denoising performance. 

5 Discussion 

The purpose of this note was to communicate the idea that one can im- 
prove (often substantially) the denoising results of NLM by replacing the 
£2 regression on patch space by the more robust £^ regression. This led 
us to propose the NLEM algorithm. The experimental results presented in 
this paper reveal two facts: (a) Phase transition phenomena - NLEM starts 
to perform better (in terms of PSNR) beyond a certain noise level, and (b) 
The bulk of the improvement comes from pixels close to sharp edges. The 
latter fact indicates that NLEM is better suited for denoising images that 
have many sharp edges. This suggests that we could get similar PSNR im- 
provements if we simply used NLEM in the vicinity of edges and the less 
expensive NLM elsewhere. Unfortunately, it is difficult to get a reliable 
edge map from the noisy image. On the other hand, observation (b) sug- 
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Table 1: Comparison of NLM and NLEM in terms of PSNR and SSIM, at noise levels a = 10, 20, ... , 100 (Averag 
over 10 noise realizations). Same parameters used: 5 = 21,^ = 7, and h = 10a. 



Image 


Method 










PSNR (dB) 










Checker 


NLM 
NLEM 


40.04 

39.73 


35.16 

34.66 


31.74 

31.66 


27.84 
29.37 


25.21 
26.71 


23.13 
24.76 


21.39 
23.22 


19.96 
21.82 


18.84 
20.50 


17.94 
19.45 


Circles 


NLM 

NLEM 


37.31 

34.27 


34.67 

34.08 


31.79 
32.33 


28.82 
30.05 


26.46 
27.92 


24.69 
26.24 


23.10 
24.87 


21.58 
23.57 


20.23 
22.36 


19.03 
21.16 


House 


NLM 


34.22 


29.78 


26.88 


25.21 


24.07 


23.37 


22.78 


22.38 


22.06 


21.81 


NLEM 


33.96 


30.10 


27.15 


25.39 


24.30 


23.51 


22.96 


22.54 


22.17 


21.95 


Barbara 


NLM 


32.37 


27.39 


24.93 


23.52 


22.64 


22.04 


21.62 


21.29 


21.07 


20.88 


NLEM 


32.11 


27.75 


25.26 


23.84 


22.90 


22.29 


21.83 


21.48 


21.20 


21.01 


Lena 


NLM 


33.24 


29.31 


27.40 


26.16 


25.24 


24.54 


24.04 


23.66 


23.34 


23.06 


NLEM 


33.15 


29.45 


27.61 


26.40 


25.53 


24.84 


24.31 


23.90 


23.53 


23.24 


SSIM (%) 


Checker 


NLM 


99.41 


98.66 


97.59 


95.51 


92.32 


88.35 


83.37 


78.05 


72.31 


66.47 


NLEM 


99.36 


98.51 


97.60 


96.30 


94.03 


91.18 


87.81 


84.04 


79.59 


74.51 


Circles 


NLM 


96.31 


94.02 


91.43 


88.66 


85.79 


82.49 


78.28 


73.28 


68.13 


63.54 


NLEM 


95.28 


93.60 


91.20 


88.74 


86.30 


83.75 


80.90 


77.85 


74.21 


70.21 


House 


NLM 
NLEM 


86.90 

86.86 


81.80 

81.25 


76.93 
77.61 


72.97 
73.85 


69.57 
70.51 


66.86 
67.77 


64.38 
65.21 


62.24 
62.96 


60.33 
60.97 


58.55 
59.03 


Barbara 


NLM 

NLEM 


89.95 
90.25 


79.48 
80.38 


71.32 
72.49 


65.20 
66.48 


60.79 
62.04 


57.42 
58.57 


54.64 
55.64 


52.58 
53.45 


50.78 
51.50 


49.27 
49.86 


Lena 


NLM 
NLEM 


86.95 
87.06 


80.41 
80.57 


76.45 
76.77 


73.42 
73.97 


70.80 
71.50 


68.45 
69.21 


66.46 
67.20 


64.66 
65.32 


62.94 
63.52 


61.32 
61.81 



gests that, by comparing the denoising results of NLM and NLEM, one can 
devise a robust method of detecting edges at large noise levels. We note 
that Wang et al. have recently proposed a non-local median-based estima- 
tor in 19]. This can be considered as a special case of NLEM corresponding 
to k = I, where single pixels (instead of patches) are used for the median 
computation. On the other hand, some authors have proposed median- 
based estimators for NLM where the noisy image is median filtered before 
computing the weights [8, 9J. In fact, most of the recent irmovations in 
NLM were concerned with better ways of computing the weights; e.g., see 
[TOllTTI . and reference therein. It would be interesting to see if the idea of 
robust Euclidean medians could be used to complement these innovations. 
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